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Abstract 

Ab initio GW calculations are a standard method for computing the spectroscopic properties 
of many materials. The most computationally expensive part in conventional implementations of 
the method is the generation and summation over the large number of empty orbitals required to 

q 

converge the electron self energy. We propose a scheme to reduce the summation over empty states 
by the use of a modified static-remainder approximation, which is simple to implement and yields 

o. ' 

accurate self energies for both bulk and molecular systems requiring a small fraction of the typical 
number of empty orbitals. 
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INTRODUCTION 



The GW methodology [l|-l3| has been successfully applied to the study of the quasiparticle 
properties of a wide range of systems {4] from traditional bulk semiconductors, insulators 
and metals to nanosystems like polymers, nanotubes and molecules The approach 

yields quantitatively accurate quasiparticle band gaps and dispersion relations from first- 
principles. A perceived drawback of the GW methodology is its computational cost; usually 
thought to be an order of magnitude more than a typical DFT calculation. One of the main 
computational bottlenecks of the traditional ab intio GW method [3j is the cost to generate 
the large number of empty orbitals needed to converge the Coulomb-hole summation term 
of the self-energy. 

Within the conventional GW approach, the quasiparticle energies and wavef unctions .{i.e., 
the one-particle excitations) are computed by solving the following Dyson equation 2|, y] (in 
atomic units): 
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where £ is the non-local, energy- dependent, self-energy operator within the GW approxi- 
mation, and E® k and ■0^ k are the quasiparticle energies and wavefunctions, respectively. In 
the typical GW approach, density functional theory (DFT) within the Kohn-Sham approx- 
imation [8j] is often chosen as the starting point for a subsequent calculation of the electron 
self-energy: the Kohn-Sham [8j wavefunctions and eigenvalues are used here as a first guess 
for their quasiparticle counterparts. In principle, Eq. [T]is a matrix equation, where E should 
be constructred in an appropriate basis. In many cases, only the diagonal elements are siz- 
able within the basis spanned by the Kohn-Sham mean-field orbitals. We assume this to be 
the case for the rest of the paper. The effects of X can thus be treated within first-order 
perturbation theory in the form £ = V xc + (S — V xc ), where V xc is the independent-particle 
mean- field approximation to the exchange-correlation potential of Kohn-Sham system . 

Within the GW approximation, the self-energy operator is expressed as £ = iGW, 
where G is the electronic Green's function and W is the dynamically screened Coulomb 
interaction. In the GW and static-COHSEX (the static limit of GW) approximations 

n n 

for the self energy, the self-energy operator, S, can be broken into two parts, [2|, |3[ 
E = S sx + Sch where Egx is the screened-exchange operator and S CH is the Coulomb-hole 
operator. The screened-exchange operator is similar to the Fock operator in Hartree-Fock 
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theory, except the bare Coulomb interaction is replaced by the screened Coulomb interac- 
tion: W GG '(q;ui) = e G Q/(q; u))v(q+ G') where v is the bare Coulomb interaction. When G 
and W are constructed in a non-self-consistent way from the DFT orbitals and eigenvalues, 
the GW approach is referred to as being within the Go^o approximation. In this article, 
all results are presented within this approximation. 

In a conventional GW calculation within the generalized plasmon-pole approximation , 
both the calculation of the Coulomb-hole self energy term: 
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and the calculation of the dielectric screening matrix, e = 1 + 47rx, at u = 0: 
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involve a summation over empty orbitals. Here N is the number of empty orbitals in the 
truncated sum, nk is a Bloch orbital with a given crystal momentum k, band index n and 
energy E n ^ v (q+ G) is the bare Coulomb interaction in reciprocal space and GG '(q) and 
^GG'(q) are plasmon-pole parameters [3]. 

There has been much research effort invested in recent years to reduce the need for empty 



orbital 
al. 



s in the GW formalism |9 



17] . The approaches by Umari et. al. [lOj and Giustino et. 



ll| eliminate the need for empty states entirely by constructing the dielectric response 



and self-e nerg y from only occupied states within a linear-response Sternheimer equation 



appraoch While these approaches eliminate the need for empty orbitals, they are 

conceptually more complicated as well as more complicated to implement and optimize. 
It is therefore still of great value to address the computational cost of the empty orbital 
generation within the traditional GW approach laid out above. 

One approach to addressing the problem of the cost associated with empty orbitals in the 
traditional GW approach is to approximate the true DFT empty orbitals with approximate 



orbitals that are computationally cheaper to generate 1 171 ]. In the recent work of 
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Samsonidze et. al. [14], the authors proposed replacing the expensive step of constructing 
the exact Kohn-Sham empty states from a traditional DFT package with a computationally 
inexpensive process of constructing the empty states from a reduced basis set consisting 
of plane- waves and resonant orbitals (generated in SIESTA 18|) orthogonalized to the real 
occupied Kohn-Sham orbitals. While it was shown that this approach vastly reduces the 
cost of generating the required empty orbitals, the approach adds significantly to the com- 
plication of the GW process. In particular, one now must run a traditional plane wave DFT 
calculation, a local orbital DFT calculation and a post-processing orthogonalization step in 
order to generate the required electron orbitals needed to proceed to the GW calculation. 
Additionally, the explicit sums in Eq. [5] and [3] must still be performed over these orbitals. 



Another approach to the empty state problem was proposed by Tiago et. al. [13j]: a 
truncation of the sum over empty orbitals in Eq. [2] can be achieved with minimal loss of 
accuracy by adding the contribution of the remaining orbitals within the static (COHSEX) 



approximation {2J, |3[. The idea relies on the fact that, in static-COHSEX, unlike GW, the 
Coulomb-hole energy can be written in a simple closed form (see next section) as well as 
in a sum over empty-states, as the static limit of Eq. [2j It was proposed that one may 
approximate the missing Coulomb-hole contribution when truncating the sum after N (i.e. 
the contribution to the sum from the empty orbitals with index between N and oo) in 
a GW calculation by their contribution to the Coulomb-hole energy in a static-COHSEX 
calculation. 



T 

al. 



lis static approximation, however, was shown to be of limited use by Bruneval et. 



12], where, instead of using the static approximation for the remaining part of the 



Coulomb-hole sum, the authors proposed using an approach based on a common non-zero 
energy denominator in Eq. [2j If a constant denominator is assumed, than one may use the 
completion relation: 

oo N 

|nk)(nk| = 1 — ^ |nk)(nk| (4) 

n=N+l n=l 

to replace the sum over the missing empty orbitals with a sum over the available orbitals. In 
order to apply this directly to Eq. [2] and Eq. [3j one must replace the n dependent denomina- 
tor with a constant. The main drawback of this common energy denominator approximation 
(CEDA) approach (known also as the extrapolar method) is that the energy denominator 
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is not uniquely denned and can only be treated as a somewhat ad-hoc parameter, and the 
_t.de energy eonvergenee is act m onotonic with this parameter. 

Recent studies by Kang and Hybertsen [19J have shown that a modified static COHSEX 
approach can be used to accurately minimize the empty orbitals problem in the Coulomb- 
hole summation of Eq [2J In that work, the authors propose completely replacing the GW 
Coulomb-hole operator with a closed-form static operator (similar to that used by Tiago 



131 ] ) with a q dependent coefficient, /(q), fit to match the GW result. This approach 
has the advantage of being completely closed form, but can be improved if one relaxes the 
requirement that the modified operator is used to replace the true GW contributions not 
only from high energy empty orbitals but also from valence and low energy conduction 
orbitals. 

In this article, we propose a modified static remainder approach based on Tiago's results 
[kJ that is more fully justified by the recent Kang-Hybertsen result [19;]. The new approach 
yields accurate GW Coulomb-hole absolute energies (to within 100 meV) with less than 10% 
of the traditionally necessary empty orbitals. Furthermore, unlike the extrapolar method 
of Bruneval |12|. this approach yields an easy to implement procedure with no adjustable 
parameters. For simplicity of presentation, we shall discuss our approach within the gener- 
alized plasmon pole model for the dielectric matrix. This approach can be straightfowardly 
applied to any existing GW computational package. 

TRADITIONAL GW CONVERGENCE WITH EMPTY ORBITALS 

As mentioned in the introduction, the band convergence of absolute energy levels in Eqh 
with the number of empty orbitals is extremely slow. In principle, one must converge a 
GW calculation with respect to the number of empty stats in both Eq. [3] and |5J However, 
for many systems, the quasiparticle energy dependence on the number of empty orbitals in 
the dielectric screening (e.g. Eq. [3]) converges much faster than with the number of empty 
orbitals in Eq. [2j For example, recent calculations for ZnO show that the Coulomb-hole 
contribution to the electronic band gap does not converge until 3,000+ empty orbitals are 
included in the summation in Eq. [2] [20[ . In Fig. (TJ we demonstrate the slow convergence of 
Eq. [2] in ZnO. The situation is even worse for nanosystems where absolute energies are often 
required for applications involving interfaces over which absolute energy level alignment is 
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needed such as the cases for molecular electronics or photovoltaic applications. 

From the bottom panel of Fig. (TJ it is immediately evident that the quasiparticle energy 
converges much more slowly with respect to the number of empty orbitals included in the 
Coulomb-hole summation, Eq. [5] than from the e summation, Eq. [3J Additionally, one 
may compute Eq. [3] in an alternative density functional perturbation theory approach that 



avoids the sum over empty orbitals. Similar techniques [ll|, l21[ to avoid the empty orbitals 
for Eq. [2] have been proposed but they are more difficult to implement and use. Therefore, 
any reduction in the summation in Eq. [2] over large numbers of empty orbitals can greatly 
reduce the cost of calculation for standard GW approaches. 



METHODOLOGY 



The static COHSEX method is the static limit of the GW approximation for the self 
energy - where everywhere e(G,G',u) is replaced by e(G, G', 0). The static remainder 
approach is based on the fact that the expectation value of the static COHSEX Coulomb- 
hole operator can be expressed either in a closed form or as a sum over empty orbitals: 



£cr>,k)= (5) 

1 N 

2 E E < nk l ei(q+G>r \ n " k ~^ ("^l e- 4(q+G,) - r ' |nk) x { [e^q; 0) - * GG ,] ^(q+G')} 

n" qGG' 



and, 

£g? /o °(n,k) = \ J2 (nk\e« G - G '> r |nk) [e^q; 0) - <W] v( q +G') (6) 

qGG' 

where N and oo denote a truncated empty state summation and closed form expression, 
respectively. Equation J5) is equal to Eq. [2] in the limit of static dielectric screening. In the 
work of Kang et. al. |l9|, the authors propose using a modified static-COHSEX operator 
that mimics the G W operator to entirely remove the need for empty orbitals. In our current 
approach, we include the full GW contribution from the low energy orbitals and add a 
single correction for the high-energy orbitals, where the static approximation is expected to 
perform well. One advantage of the present approach is that it can be used in conjunction 




FIG. 1: Top: The convergence of the Coulomb-hole contribution to the self-energy, Eq. [21 with 
respect to the number of orbitals included in the summation, N, using a dielectric matrix calculated 
with 1000 empty bands. For all calculations on ZnO, a 5x5x4 k-point grid is used. Bottom: The 
convergence of the quasiparticle energy, Eqp, with respect to empty states in the polarizability 
sum Eq. [3] and with respect to empty states in the Coulomb-hole sum Eq. [21 The red curve shows 
the VBM Eqp in ZnO using a fixed 3,000 bands in the Coulomb-hole summation and varying the 
number of bands included in the polarizability summation. The black curve shows the VBM Eqp 
in ZnO using a fixed 1,000 bands in the polarizability summation and varying the number of bands 
included in the Coulomb-hole summation. 

with a full-frequency (as opposed to GPP model) screening approach to both calculate 
the fine structure of energy dependence of the self energy, E(tu), as well as converging the 
absolute value with respect to empty orbitals. In our modified static remainder approach, 
we calculate both the GW Sch partial sum (Eq. [2]) and the COHSEX Ech partial sum (Eq. 
[5]) up to the number of DFT bands available. We then add a modified static correction to 
the GW Coulomb-hole energies: 
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(nk|E£ H (r,r';£)|n'k)= (7) 
(nk| ££ H (r, r'; E) \n'k) + V (nk| S^ /oo (r, r') |n'k) - (nk| S^ /JV (r, r') |n'k> ) . 

The factor of 1/2 in Eq. [7] is justified from the recent work of Kang and Hybertsen 19), 
where the authors show that the GW contribution of high energy bands (corresponding to 
large G- vectors) to the Coulomb-hole self energy asymptotes to 1/2 of the equivalent static 
COHSEX band contribution. 

One may qualitatively derive this result from Eq. [21 if one assumes that we are interested 
in a state n with energy, E, near zero, and that, for a given high n", the sum over matrix 
elements are dominated by a small set of q + G near such that |^(q+ G)| 2 /2m « E n » 
and that the plasma frequency for large q and G obeys a homogeneous gas dispersion 
^g,g (q) ~ |^(q + G)| 2 /2m. In this case, the contribution to Eq. |2] from a high-energy 
empty orbital n" reduces to: 



(nk| ££ H (r, r'; E) \n'k) = - ( nk l e i(q+G) ' r |n"k-q) (n"k-q| e -'^+ G ')- r ' | n 'k) (8) 
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j ( nk l e i(q+G) ' r |n"k-q) (n"k-q| e - i(q+G,) - r ' \nk) 
x{[eGG'(q;0)-W] u(q+G')}, 
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which, when compared to Eq. \5[ confirms the factor of ~. 

The Coulomb-hole self-energy contribution to the convergence of energy levels in bulk 
silicon (using a 5x5x5 k-point grid) and the silane (SiH 4 ) molecule (in a supercell calculation) 
are shown in Fig. [2] as a function of the band cutoff, N, in Eq. [2] The convergence on energy 



levels in silane is significantly slower than that in silicon [22] because of the large number 
of free-electron like vacuum states. The silicon calculations were done with a 25 Rydberg 
wavefunction cutoff and a 10 Rydberg dielectric matrix cutoff. The silane calculations were 
done with a 75 Rydberg wavefunction cutoff and a 6 Rydberg dielectric matrix cutoff. The 
needed volume of the supercell used, (25cm) 3 , and the corresponding number of vacuum 



states, is minimized by using a truncated Coulomb interaction 



23( . Despite this, the largest 




# of Bands in Partial Sum # of Orbitals in Partial Sum 

FIG. 2: Comparison between the contributions to the Coulomb-hole sum for the full GW operator 
vs. results from the 1/2 the static COHSEX Coulomb-hole operator for orbitals beyond the number 
of real DFT bands/orbitals used: 12 in silicon and 100 in Silane. A 5x5x5 k-point grid is used in 
Si. The plotted quantity is J2n"=n DFT +i E q GG' ( nk l e i( .^+ G > r |n"k-q) (n"k-q| e - l ^+ G ')- r ' |n'k) x 
Igg,(q, n, n' , n") where I CH is the term in {} in Eqs. ([2]) and ([5]) respectively. 

computational cost in the GW calculation on silane is the DFT generation of the empty 
orbitals, representing more than 50% of the total computational expense. The calculation 
of the polarizability and the evaluation of the self energy require less computational time, 
and they scale nearly linearly to thousands of CPUs. 

RESULTS 

A comparison between the convergence of the residual value of the GW expression (Eq. 
[2]) and 1/2 of the static COHSEX approximation (Eq. [5]) for the Coulomb-hole contribution 
to the electron self-energy starting at some u-dft is shown in Fig. [2] for the valence band 
maximum (VBM) in Silicon and the highest occupied molecular orbital (HOMO) in the 
silane molecule. The figure shows the cumulative contributions of the high-energy orbitals 
to Sqh for both the GW operator and the 1/2 static COHSEX operator for orbitals above 
12 and 100 for silicon and silane, respectively. The residual value of the 1/2 static COHSEX 
results reproduce the equivalent GW curves extremely well. Therefore, replacing the GW 
operator with the modified static remainder in Eq. [7] yields very good agreement with a 
fully converged GW calculation. This justifies the truncation of the partial sum in Eq. [2] 
and the addition of the modified static remainder correction. For both silicon and silane, 
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FIG. 3: Coulomb-hole energies of the valence band maximum in Si (left) and ZnO (right) in 
the modified static-remainder approach compared to the energies from the standard approach of 
truncating the Coulomb-hole summation in Eq. [2] as a function of the number of DFT bands. In 
the static-remainder approach the summation is also truncated at the same number of bands but 
the modified static remainder is added to the sum. A 5x5x5 and 5x5x4 k-point grid is used in Si 
and ZnO respectively. The grey lines represent the result using the maximum number of bands 
and the static remainder included. 

one can get a converged £ch to within 100 meV with less than 10% of the original number 
of empty orbitals required. This is a very high level of accuracy considering the modified 
static correction in both cases is greater than 1 eV. Even higher accuracy may be reached if 
one increases the number of actual DFT empty orbitals used. In the case of Si, an accurate 
Eqh (within 100 meV) can be reached with the use of only 10 empty bands when a 5x5x5 
k-point grid was used. Furthermore, the convergence for this approach is nearly monotonic 
in terms of the number of empty Kohn-Sham orbitals employed in the calculation. Figure [3] 
shows the convergence of the modified static remainder corrected Sch and the uncorrected 
GW S CH as a function of DFT empty bands used for Si an ZnO - the grey line corresponds 
the "best guess" final value using the static-remainder on top of the largest number of real 
DFT bands. 

In table I, we show the convergence behavior with respect to empty orbitals of our 
GW+static-remainder approach compared to a traditional GW approach for bulk Si, MgO, 
ZnO, and solid Ar. In all cases, the modified static- remainder approach significantly im- 
proves the convergence rates. An accuracy of less than 100 meV in the absolute energies 
can be typically reached with only a few conduction bands. 
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Si - r 


10 Bands 


40 Bands 


80 Bands 


160 Bands 


n=l 


-5.99 


6.30 


6.47 


6.53 


n=l (SR) 


-6.56 


6.55 


6.55 


6.55 


n=4 


6.19 


5.62 


5.30 


5.15 


n=4 (SR) 


4.99 


4.99 


5.04 


5.08 


n=5 


9.46 


8.90 


8.60 


8.48 


n=5 (SR) 


8.35 


8.33 


8.39 


8.41 


n=10 


15.23 


14.48 


14.10 


13.94 


n=10 (SR) 


13.76 


13.73 


13.80 


13.84 


ZnO - r 


100 Bands 


500 Bands 


1500 Bands 


3000 Bands 


n=26 


6.11 


4.48 


4.00 


3.90 


n=26 (SR) 


3.88 


3.73 


3.80 


3.82 


n=27 


7.97 


7.32 


7.22 


7.21 


n=27 (SR) 


7.22 


7.19 


7.20 


7.21 


MgO - r 


50 Bands 


200 Bands 


450 Bands 


900 Bands 


n=4 


-2.09 


-2.86 


-2.95 


-2.96 


n=4 (SR) 


-3.15 


-3.02 


-2.97 


-2.97 


n=5 


5.10 


4.81 


4.78 


4.78 


n=5 (SR) 


4.70 


4.77 


4.78 


4.78 


Ar - r 


50 Bands 


150 Bands 


375 Bands 


750 Bands 


n=4 


-7.64 


-8.19 


-8.39 


-8.42 


n=4 (SR) 


-8.52 


-8.50 


-8.43 


-8.43 


n=5 


5.26 


5.38 


5.411 


5.42 


n=5 (SR) 


5.45 


5.43 


5.42 


5.42 



TABLE I: Convergence of the Eqp (in eV) within the GqWq approximation with respect to the 
number of DFT orbitals used in the Coulomb-hole summation for several material systems. Here 
n refers to the band index and (SR) refers to the addition of the static remainder. 
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FIG. 4: Coulomb-hole part of the self-energy, with and without the static-remainder, for the highest 
occupied molecular orbital (HOMO) of the BND (bithiophene naphthalene diimide - shown in inset) 
molecule as a function of the number of DFT orbitals included in the Coulomb-hole sum. The grey 
line represents the result using the maximum number of bands and the static remainder included. 



To test the modified static reminder approach on a large molecular system, we compute 
the Coulomb-hole contribution to the self-energy for the BND (bithiophene naphthalene 
diimide) molecule containing 46 atoms 24J. The supercell was set to 76.93 x 36.31 x 20.18 
atomic units. The calculations were done with a 60 Rydberg wavefunction cutoff and a 
6 Rydberg dielectric matrix cutoff. The polarizability was computed with 953 orbitals 
(78 occupied + 875 empty orbitals up to 1 Rydberg cutoff in DFT eigenvalues), and the 
Coulomb-hole part of the self-energy was evaluated as a function of the number of orbitals, 
as shown in Fig. HI One can see that the Coulomb-hole term computed with 953 orbitals 
without the addition of the remainder is only converged to within 1 eV. Including the static 
remainder correction improves the convergence to better than 0.1 eV. 

In conclusion, we have presented a modified static remainder approach that reduces 
the number of empty states involved in evaluating £ch by over an order of magnitude. 
This approach is particularly useful when applying the GW method to molecules and other 
nanostructures where absolute energies, as opposed to just energy gaps, are desired. A 
limitation of this method is that it does not address the problem of the sum over empty 
states required in evaluating the dielectric matrix (Eq. [3]) . However, the dielectric matrix 
converges faster than the absolute energies of Ech for many solids 20[ and can more easily 



be replaced by calculation using the density functional perturbation theory approaches. Our 
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approach here shows nearly monotonic convergence towards the converged GW Sch values 
and can be implemented in a simple and automatic way in standard GW computer codes. 
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